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ABSTRACT 

We present a formulation for multigroup radiation hydrodynamics that is correct to order 0(y/c) 
using the comoving- frame approach and the flux-limited diffusion approximation. We describe a 
numerical algorithm for solving the system, implemented in the compressible astrophysics code, CAS- 
TRO. CASTRO uses an Eulerian grid with block-structured adaptive mesh refinement based on a 
£SJ ■ nested hierarchy of logically-rectangular variable-sized grids with simultaneous refinement in both 

space and time. In our multigroup radiation solver, the system is split into three parts, one part that 
^ . couples the radiation and fluid in a hyperbolic subsystem, another part that advects the radiation in 

frequency space, and a parabolic part that evolves radiation diffusion and source-sink terms. The hy- 
" perbolic subsystem and the frequency space advection are solved explicitly with high-order Godunov 

schemes, whereas the parabolic part is solved implicitly with a first-order backward Eulcr method. 
Our multigroup radiation solver works for both neutrino and photon radiation. 
Subject headings: diffusion - hydrodynamics - methods: numerical - radiative transfer 

1. INTRODUCTION 

fS j Numerical simulations arc a useful tool for many radiation hydrodynamic problems in astrophysics. However, 
numerical modeling of radiation hydrodynamic phenomena is very challenging for a number of reasons. To achieve a 
highly accurate description of the system usually requires a multidimensional treatment with high spatial resolution. 
It is also desirable to handle the radiation frequency variable properly, because radiation transport is a frequency- 
^ 1 dependent process. Moreover, the numerical algorithms for radiat ion hydrodynamics m ust be stable and efficient. 

Co re-collapse supernovae are such complex phenomena (jBethei 119901 : iKotake et al.l 120061 : iJanka et al.l l2007t Uankal 
l2012f ) that neutrino radiation hydrodynamics simulations have been indispensable for helping understand the explosion 
mechanism. Various algorithms have been developed to solve neutrino radiation hydrodynamics equations. Because 
^ | of the complexities involved in the problem and the limited computer resources available , various tradeoffs have 
I/"*) . to be made. Many of these algorithms are ID (e . g.. iBowers fc Wilson! ITpSl 1 ruennl [19851 : iMezzacappa fe Bruennl 
^" ■ 119931: iBurrows et all 120001: iLiebendorfer et al.l 12004 lO'Connor fe Ottj|2010Tr"or have used the so-called "ray-by-ray" 
\ approach, which does not truly perform multi-dimensional transport (e.g. Burrows et al.l lT995: Ra mpp fe Jankall2002l : 
. iBuras et al.ll2006l : lTakiwaki et al.ll2012D . The two-dimensional code of lMiller et al.l ( 19931) is a gray flux-limited diffusion 



o 



£ — , ■ (FLD) code that does not incl ude order Q ( v/c) t erms, where v is the characteristic velocity of the system. The 2D 
' multigroup multiangle code of iLivne et al.l (|2004l ) does not include all order Q(v/c) terms, and in parti cular it lacks 
£N| . direct coupling among radiation groups. The two-dimensional algorithm of iHubenv fc Burrows! (|2007l ) is based on 
a mixed-frame formulation of a mu ltigroup two-moment system that is correct to order 0(v/c), but it remains to 
be implemented in a working code. iSwestv fc Mvral (|2009f ) have developed a fully coupled 2D multigroup FLD code 
based on a comoving frame formulation that is correct to order 0(v/c). Besides the grid-based methods, smoothed 



particle hydrodynamics (SPH) methods have also been applied to multi-dimensional neutrino radiation h ydrodynamics 
5^ 1 simulations of core-collapse supernovae (e.g.. lHerant et al.l 1994HFrver fc Warren! [gOOllFrver et~alll2006l . but the SPH 



& , scheme has not yet been extended to multigroup. Recently 



Abdikamalov et al.1 (|2012l ) developed a Monte Carlo method 



for neutrino transport, but hydrodynamics and radiation transport are not yet coupled in the scheme. 

In this paper, we present a numerical algorithm for solving multi-dimensional multigroup photon and neutrino 
radiation hydrodynamics equations in the comoving frames that are correct to order 0{v/c). Radiation quantities 
that are of interest are the radiation energy density, radiation flux, and radiation pressure tensor. It is more accurate 
to evolve both the radiation energy density and the radiation flux especially when the radiation is highly anisotropic 
and optically thin. However, the two-moment approach is very expensive in terms of computer time and memory for 
multi-dim ensional multigroup ra diation hydrodynamics simulations. Thus, we have adopted the flux-limited diffusion 
approach (jAlme fc Wilson l ll9 73j for computational efficiency. In the FLD approach, only the radiation energy density 
needs to be evolved over time, and the radiation flux is derived through a diffusion approximation. Furthermore, we 
use an analytic closure for the radiation pressure tensor. We have adopted a multigroup metho d, in which the radiatio n 
frequency is discretized into multiple groups. The equations we solve are similar to those in ISwestv fc Mvral <|2009h . 
However, we have neglected inelastic scattering of neutrinos on matter, which is currently treated as elastic scattering. 
Since the total cross section for the latter at the neutrino energies near and outside the neutrinospheres are ~ 20-100 
times smaller than the cross sections for scattering off nucleons and nuclei, not including inelastic effects at this stage 
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has a minor effec t on th e qualitative behavior of collapse, bounce, and shock dynamics. Nevertheless, the approach of 
iThompson et al.l (|2003f) to handle inelasticity an d energy redistribution can be incorporated into our scheme. Another 
difference between our scheme and the scheme of iSwestv fc Mvral (I2009D is that we have developed a Godunov scheme 
in which radiation is fully coupled into a characteristic-based Riemann solver for the hyperbolic part of the system. 

We have implemented our algorithm in the compressible astrophysics code, CASTRO. This is the third paper in 
a series on the CASTRO code and its numerical algori thms. In our previou s papers, we describe our treatment of 
hydrodynamics, including gravity and nuclear reactions (|Almgren et al.ll2010l hen ceforth Paper I), a nd our algorithm 
for flux-limited gray radiation hydrodynamics based on a mixed-frame formulation (jZhang et al.l20lTl henceforth Paper 
II). Here, we describe our algorithm for flux-limited multigroup radiation hydrodynamics based on a comoving- frame 
formulation. 

CASTRO uses an Eulerian grid with block-structured adaptive mesh refinement (AMR). The AMR algorithm in 
CASTRO uses a nested hierarchy of logically-rectangular variable-sized grids and refines simultaneously in both space 
and time. A major advantage of CASTRO is that it can efficiently achieve high resolution through the use of AMR. 

In § [2l we present the multigroup radiation hydrodynamics equations that CASTRO solves. We show that the 
system can be divided into a hyperbolic subsystem (§ 12. 3|) . a subsystem of advection in frequency space (§ I2.4|) . and 
a parabolic subsystem for radiation diffusion (§ 12. 5|) . Analytic results for the mathematical characteristics of the 
hyperbolic subsystem are also presented in § 12.31 The hyperbolic and the frequency space advection steps are solved 
by an explicit solver, whereas the parabolic subsystems is solved by an iterative implicit solver. We describe the 
explicit solvers in § 13.11 followed by a description of the implicit solver in § 13.21 We also discuss acceleration schemes 
that can greatly improve the convergence rate of the iterative implicit solver (§ I3.2.4[) . In §[4j we present results from 
a series of test problems. A scaling test is shown in § [5] to demonstrate the scaling behavior of the multigroup group 
radiation solver. Finally, the results of the paper are summarized in § [6] 

2. MULTIGROUP RADIATION HYDRODYNAMICS 

In this section, we present the multigroup radiation hydrodynamics equations that CASTRO solves. Our formulation 
is based on a comoving-frame approach, unlike our gray radiation hydrodynamics solver presented in Paper II. In the 
comoving-framc approach, the radiation quantities and the opacities are measured in the comoving frames. The set 
of equations are correct to order 0(v/c). We then derive the multigroup radiation hydrodynamics equations using 
frequency space discretization. We also analyze the mathematical characteristics of the system. Here we present the 
results in terms of neutrino radiation and indicate modifications needed for photon radiation. 

2.1. Frequency- Dependent Radiation Hydrodynamics Equations 
The system of neutrino radiation hydrodynamics that is correct to order 0(v/c) can be described by 

dp 



dt 



(pu) = 0, (1) 



di " PU) 1 V • (puu) + Vp = r^F v dv + F G , (2) 



dt J c 

+ V ■{pEu+pu)=u- I ^F, 
d{ P Y e 



V ■ (pY e u) = / £{ckE v - 47T V )diy, (4) 



dt 

^ + V • (E v u) = - V ■ F v - {ckE v - Airn) + g^(P v ■ Vu), (5) 

(|Buchlerlll98"llMihalas fc Mihalasl[T999t lCastoril2004l : ISwestv fe Mvrdl2009| ) . Here p, it, p, and E are the mass density, 
velocity, pressure, and total matter energy per unit mass (internal energy, e, plus kinetic energy, u 2 /2), respectively. 
Fq is the external force on the fluid (e.g., gravitational force). E u , F U1 and P„ are the monochromatic radiation 
energy density, radiation flux, and radiation pressure tensor at frequency v, measured in the comoving fluid frame, 
respectively. The speed of light is denoted by c. The absorption coefficient is denoted by k, and the total interaction 
coefficient including conservative scattering is given by x- The emission coefficient is denoted by ry. For photon 
radiation £ is zero, whereas for neutrinos £ is given by 

hv 

where me is the baryon mass, h is the Planck constant, and s = 1 for v e neutrinos, s = — 1 for v e neutrinos, and s = 
for other flavor neutrinos. Also note that for neutrino radiation with multiple species, the integrals over frequency 
imply summation over the species. The colon product : in P„ : Vit indicates contraction over two indices as follows, 

P v : Vu = (P v ).. (Vu) y . (7) 

The system of the radiation hydrodynamics (Equations Q] [5]) is closed by an equation of state (EOS) for matter, the 
diffusion approximation for radiation flux, and a closure relation for radiation that computes the radiation pressure 
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tensor from radiation energy density and radiation flux. In the FLD approximation (lAlme fc WilsonlH97l . the 
monochromatic radiation flux is written in the form of Fick's law of diffusion, 

F V = -—VE V , (8) 

X 

where A is the flux limiter at frequency v. We have implement ed a variety of choices for the flux limiter (e.g., 
Minerbol fl978t IBruenn etHI fl~978t E evermore fc Pomranind [1981|) . For example, one form of the flux limiter in 



Levermore fc PomranineT i 19810 is 

A = 2 + R (9) 
6 + 3R + R 2 ' () 

where 

In the optically thick limit, R — > 0, A — > 1/3, and the flux approaches the classical Eddington approximation F v « 
— (c/3x)V£' l /. In the optically thin limit, R — > oo, Aw 1/iJ, and F v ft '■ cE v as expected for s treaming radiation. For 
the radiation pressure tensor, we use an approximate closure (jMincrbolll978t lLevermo"relll984D 

P» = \0- - f)EJ + ^(3/ - l)E v nh, (11) 

where 

/ = A + X 2 R 2 (12) 

is the scalar Eddington factor, I is the identity tensor of rank 2, and n = V-E„/|V-E„|. There are other choices for the 
scalar Eddington factor. For example, iKershawl (|1976|) has suggested 

f = \ + \x 2 R 2 . (13) 

Note that in the optically thick limit, / — > 1/3, and in the optically thin limit, / — > 1. 

Using the FLD approximation and the approximate closure, the equations of the radiation hydrodynamics now 
become 

^+V-(pu) = 0, (14) 
1 V • {puu) + Vp= - ( \VE v dv + F G , (15) 



dt 
d( P E) 



, V • (pEu +pu)= — u • / XVE u dv + / (ckE u - Airiri&v + u- F G , 
dt Jo Jo 

+ V • (pY e u) = j t(cKE u - 

BE / cX \ 

+ V • {E u u) = V • f — VE„ J - {ckE u - 



(16) 
(17) 



1 -J-E V XI -u+ 3 f 1 E v nh :Vm). (18) 



d In v 



2.2. Multigroup Radiation Hydrodynamics Equations 

To numerically solve Equations (fT4")) - (fT5)) . we first discretize the system in frequency space by dividing the radiation 
into N energy groups. We define the radiation energy density in each group as 

E g = / E u du, (19) 

where P g -x/2 and P g +i/2 are the frequency at the lower and upper boundaries, respectively, for the g-th group. For 

clarity, J" 9+ ^ d^ is denoted by J g Au hereafter. 

In the multigroup methods, various group mean coefficients are introduced to the system. For example, the Rosscland 
mean interaction coefficient is defined as 

_ J gX -HdB„/dT)4» 
Xg f(dB v /dT)dv ' { ' 
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where B v = tj/k. Unfortunately, it can be costly to compute these group mean coefficients because the interaction 
coefficients usually depend on frequency However, the difference between various group means will diminish as more 
groups are used. Hence we will simply assume that 



Xg = X{Vg)- 



(21) 
(22) 



where v g is the representative frequency for the 5-th group. We also assume that the flux limiter is approximately 
independent of frequency within each group, and denote the flux limiter and the scalar Eddington factor for the p-th 
group as X g and f g , respectively. As for the emissivity, we define 



Jg 



47T 

—ri(i>g)Ai/ g , 



(23) 



where Ai/ g = v g +\/2 — Vg-1/2 is the group width. However, in the case of photon radiation, we have adopted the 
polylogarithm function based method of iClarkl (|1987t ) for computing the group-integrated Planck function. Thus, 
assuming local thermodynamic equilibrium, we have for photon radiation 



Jg 



Air 
c 



where B g = J g B v Av. 

Using these definitions and approximations, we obtain the multigroup radiation hydrodynamics equations 

dp 
at 

d(pu) 



V • (pu) = 0, 

, V • (puu) + Vp + J2 X 9^E g = F G 

g 

g 

d( P Y e ) 



Of 



dEg 

dt 



3-/ s 



E„u 



dt 

u ■ V 



Y^C(KgEg -jg)+U- F G , 

g g 

V-( P Y e U)= Y,C£g(KgEg 



Jg) 



l-/ 9 



E, 



c(K g Eg ~ jg) + V 



C\g 

Xg 



■V£ £ 



(24) 

(25) 
(26) 

(27) 

(28) 

(29) 



0_ 

dv 



1 — f 3 f — 1 

^-V • u + — nn : Vti I vE v 



Av — E n nn : Vu, 

n y 



where g = 1, 2, . . . , N, £ g = £(y g ), J2 q denotes Y^=i 



2.3. Mathematical Characteristics of the Hyperbolic Subsystem 

The system of multigroup radiation hydrodynamics (Equations 1251 129[) . similar to the system of gray radiation 
hydrodynamics we have studied in Paper II, has a hyperbolic subsystem, 

dp 
dt 

d(pu) 



dEg 

dt 



V ■ {pu) = 0, 

, V • (pun) + V P + Y, A S V£ 9 = F G , 
g 

V • (pEu + pu) + u ■ \VE g = u -F G , 

g 

d(pY e ) 

'3-/„ 



(30) 
(31) 



dt 



dt 



V • 



■E g u I — u ■ V 



V • (pY e u 
1-/, 



E a = 



(32) 

(33) 
(34) 



where g = 1, 2, . . . , N. Since Y e is just an auxiliary variable that does not alter the characteristic wave speeds of the 
hyperbolic subsystem, we will neglect it in the remainder of this subsection for simplicity. Gravitational force terms 
are source terms, and thus they are not included in the following analysis. The term — [(3/ g — 1) /2]E g nh : Vu in 
Equation (|29|) is not included in the analysis of the hyperbolic subsystem here even though it becomes comparable to 
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the term V- {[(3 — f g )/2}E g u} in Equation (fM)) in the streaming limit where f g —tl. The first reason for not including 
the term is that the main purpose of including radiation in the hyperbolic subsystem is to improve the accuracy in the 
optically thick limit where the radiation can affect the hydrodynamics significantly. The term — [(3/ ff — l)/2]E g hii : Vit 
is not significant in the optically thick limit where f g — > f/3. The second reason is that in the streaming limit, where 
A — > 0, radiation force has little impact on the motion of matter. Thus, radiation is essentially decoupled from the 
hydrodynamics of matter in the streaming limit. Moreover, in the streaming limit the diffusion term V • [(cA 9 /x s )V E g ] 
(Equation [25]) will dominate the term — [(3/ g — l)/2]E g hn : Vu and all the terms in Equation Therefore, the 

term — [(3/ g — l)/2]E g hh : Vit is insignificant in both the diffusion limit and the streaming limit. 

In CASTRO, we solve the hyperbolic subsystem with a Godunov method, which utilizes a characteristic-based 
Riemann solver. The Godunov method requires that we analyze the mathematical characteristics of the hyperbolic 
subsystem. For simplicity, let us consider the system in one dimension, which can be written in terms of primitive 
variables as, 

dQ 
dt 



A 8Q 

OX 



0. 



(35) 



where the primitive variables are 



Q = 



( p \ 

u 

p 

Ei 
E 2 



\E N j 



(36) 



and the Jacobian matrix is 



.4 



(u 
' 






p 

u 

IP 

£i(3-/i)/2 
£ 2 (3-/ 2 )/2 





1/P 
u 







VO E N (3-f N )/2 





Ai/p 

u 










A 2 /p 









\ 

Xn/p 






(37) 



where 7 is the adiabatic index of the matter. This system is hyperbolic because the Jacobian matrix is diagonalizable 
with N + 3 real eigenvalues, 

u — c s , u + c s , u, u, u, . . . , u. (38) 



Here 



\ 



P 

9 



3-/ fl 



\ E 9 




(39) 



is the radiation modified sound speed, where c m = \/yp/p is the sound speed without radiation and c g = 

\/[(3 /s)/2](A 3 £ , g/p) is the contribution to the sound speed from each group. The corresponding right eigenvectors 
are, 



-C-slP 
1P/P 

cl/A 2 



( / ^ 

Cs/P 

ipIp 
4/x 1 

4/X2 



V c 2 N /X N J V c n/Xn J 
and the corresponding left eigenvectors are, 



(0. 
(0. 

(1. 
(0. 
(0. 

(0. 



-p/2c s 
p/2c s 



0. 
0. 
0. 



f/2c, 2 , 
l/2c s 2 , 

-wi/Ai, 
-W2/A2, 









( ° \ 





-Ai 
1 




Ai/2c. s 2 , 
Ai/2c s 2 , 
-A!/c s 2 , 

1 — W\, 

-W2A1/A2, 



/ ° \ 



V J 





-A 2 

1 



A 2 /2c s 2 , 
A 2 /2c s 2 , 
-A 2 /c s 2 , 
-W1A2/A1, 

f - ^2, 



0, —ujn/Xn, — wjvAi/Ajv, — watA 2 /Ajv, 



V J 



( ° \ 






—Xn 





V 1 / 



Aat/2c s 2 ) 
A w /2c s 2 ) 
-X N /c s 2 ) 
-wiAiv/Ai) 
-w 2 Ajv/A 2 ) 



(40) 



(41) 
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where 



I- (42) 



These N + 3 eigenvectors define the characteristic fields for the one-dimension al system. By computing the product 
of the right eigenvectors and the gradients of their corresponding eigenvalues (|LeVequel l2002f ) . we find that the first 
and second fields are genuinely nonlinear corresponding to either a shock wave or a rarefaction wave. The rest of the 
fields are linearly degenerate corresponding to a contact discontinuity in either density or radiation energy density 
of one group moving at a speed of u. Note that there is also a jump in fluid pressure such that the total pressure, 
Ptot — P + Sg-^s^t/; i s constant across the contact discontinuity. Obviously, in three dimensions, there are two 
additional contact discontinuities for transverse velocities just like the case of pure hydrodynamics. 

2.4. Advection in Frequency Space 
The integrals in Equation (|!?9")) form the following system, 



8Eg 

dt 



d_ 
dv 



3/-1 



fin : Vu uE, 



dv, 



which can also be written in the conservation law form 



dt 



d{a q qi) 
dim; 



= 0, 



(43) 



(44) 



where qi = vE v and a q = — [(1 — /)/2](V • u) — [(3/ — l)/2]nn : Vit. Here we have used Equation (fl"9"]) . Equation (|44|) 
describes the advection of q\ in logarithmic frequency space with an advection speed of a q . Because / gidlni/ is the 
integrated radiation energy, the system conserves energy. 

Various terms in Equation (|29[) have been split between the hyperbolic subsystem (§ I2.3[) and the system for the 
advection in the frequency space (§ 12 .4[) . There is however another way to look at Equation (|29l) . The evolution of 
the radiation energy density (Equation 129 j) without the diffusion term V • [(c\g/xg)^E g ] and the source-sink term 
—c(K g Eg — jg) can also be split as 



dE^ 

dt 
dE v 

dt 



V ■ (EgU), 



8 



din v 



3/-l„ 



fin : Vm E v 



2 2 

Equation (|46[) also represents an advection in logarithmic frequency space and can be written in the form 

dq 2 d{a q q 2 ) 



dt 



dlnu 



0. 



(45) 
(46) 

(47) 



where q 2 = E v . Here the conserved quantity is the number of radiation particles (i.e., J E v jvAv). This new way of 
splitting will be further discussed later in Section l3~Tl 



2.5. Multigroup Radiation Diffusion 

The parabolic part of the system consists of the radiation diffusion and source-sink terms, which were omitted from 
the discussion of the hyperbolic subsystem (§ I2.3|) and advection in frequency space (§ 12.4)1 . 

(48) 
(49) 

(£v4 (50) 

where e is the specific internal energy of the fluid, and g = 1, 2, . . . , N. The term c{n g Eg — j g ) represents the energy 
exchange in the comoving frame between the matter and the g-th radiation group through absorption and emission of 
radiation. Note that the terms on the right-hand side (RHS) of Equations (|48p. (|4"9")l . and (|5U)) all contain the speed of 
light c. Thus these are order O(l) terms in v/c. An implicit treatment for these equations is usually necessary because 
of their stiffness. For photon radiation, Equation (|49|) is not needed. 



d(H ^ 

9 

= 2-j ^aKKgEg -jg), 



m 

dEg 

dt 



= -c(KgEg - jg) + V 
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3. NUMERICAL METHODS 

CASTRO uses a nested hierarchy of logically-rectangular, variable-sized grids with simultaneous refinement in both 
space and time. The AMR algorithms for hyperbolic and radiation diffusion equations in CASTRO have been described 
in detail in Papers I and II, respectively. The multigroup solver does not require any changes in the AMR algorithms. 
Thus, we will not describe them again in this paper. 

In this section, we describe the single- level integration scheme. For each step at a single level of refinement, the 
state is first evolved using an explicit method for the hyperbolic subsystem and advection in frequency space (§ 13. 
Then an implicit update for radiation diffusion and source-sink terms is performed (§ 13. 2[) . In § 13.2.41 we describe 
two acceleration schemes for speeding up the convergence rate of the implicit solver for multigroup radiation diffusion 
equations. 



3.1. Explicit Solver for Hyperbolic subsystem and Advection in Frequency Space 

The hyperbolic subsystem (§ 12. 3[) is treated explicitly. Our Godunov algorithm for the hyperbolic subsystem of the 
multigroup radiation hydrodynamics is an extension of the hyperbolic solver for the gray radiation hydrodynamics 
presented in Paper II , which in turn is based on the hydrodynamics algorithm presented in Paper I of this scries. 
We refer the reader to Papers I and II for a detailed description of the integration scheme, which supports a general 
equation of state, self-gravity, nuclear reactions, the advection of nuclear species and electron fraction. Here we will 
only present a brief overview of the Godunov scheme. 

The time evolution of the hyperbolic subsystem can be written in the form 

f-V.F, (5!) 

where U = (p, pu, pE, E±, E2, . ■ . , En) t denotes the conserved variables, and F is their flux. The conserved variables 
arc defined at cell centers. We predict the primitive variables, including p, m, p, pe, E g , where g — 1,2,..., A, 

from cell centers at time t n to edges at time t n+1 / 2 and use an approximate Riemann solver to construct fluxes, 
F ™+ 1 / 2 j on cell faces. This algorithm is formally second order in both space and time. The time step is computed 
using the standard Courant-Fricdrichs-Lewy (CFL) condition for explicit methods, and the sound speed used in the 
computation is now the radiation modified sound speed c s (Equation 139 j). The time step is also limited by maximally 
allowed changes in temperature and electron fraction in the implicit solver. Additional constraints for the time step 
are applied if additional physics, such as burning, is included. CASTRO solves the hyperbolic subsystem of radiation 
hydrodynamics with an unsplit pieccwisc parabolic method (PPM) with characteristic tracing a nd full cor n er cou pling 
dMiller fc Colella] |2002f) . An approximate Riemann solver based on the Riemann solver of iBell et al~l (|1989f ) and 
IColella et all (|1997l T that utilizes the eigenvectors of the system, is used to compute the Godunov states and fluxes at 
the cell interface. 

Advection in frequency space (§ 12. 4[) can also be treated explicitly because its time scale is comparable to that 
of the hyperbolic subsystem. One approach is to adopt an operator-splitting meth od and to evolve Equation (|44"|) 
with a Godunov method. This is the approach taken by Ivan der Hoist etaLl (|2011l ) in the CRASH code for photon 
radiation. The term — [(3/ ff — l)/2\E v iifi : Vit can be treated as a sou rce term of Equation (l44l. The success of 
this approach for a number of test problems has been demonstrated bv Ivan der Hoist et al.l (|2011l ). We have also 
found that this approach can produce satisfactory results for all the photon test problems in Section [4] However, 
we have found that erroneous results were obtained for core-collapse supernova simulations (§ 14.21) . probably due to 
operator-splitting errors. One difficulty in core-collapse supernova simulations is that the divergence of velocity can be 
significant because there is a strong shock and the typical velocity in front of the shock is ~ 0.1-0.2 c. In the case of 
core-collapse supcrnovae, the negative divergence of velocity will shift neutrinos towards high frequency (Equation [44]). 
The stronger conservative scattering at higher frequency can further amplify the neutrino energy density due to the 
work term u ■ V{[(1 — f g )/2]E g } fEquation 134)) . These two processes can combine to cause unphysical results due to 
splitting errors. However, the origin of the various terms involved in these two processes is the term d(P u : Vit) / 'dhiis 
in Equation (|5|). To avoid operator-splitting errors, it is therefore desirable to treat the term d(P v : Vu) / d\nv without 
splitting it. On the other hand, the term can also make significant contribution to the hyperbolic subsystem in the 
optically thick limit, and we do not wish to sacrifice accuracy of the hyperbolic subsystem by taking the traditional 
approach of leaving out radiation from the Riemann solver for the hydrodynamics. 

To avoid the operator-splitting without sacrificing the accuracy of the hyperbolic subsystem, we have adopted a 
new approach based on splitting the radiation energy density equation into Equations (|45l) and (|46|) . Note that the 
radiation diffusion and source-sink terms will treated separately by an implicit solver (§ 13. 2[) . Our new approach 
still utilizes the Godunov scheme for the hyperbolic subsystem governed by Equations (|3tj)) - (|34p . After the Godunov 
states and fluxes at the cell interfaces are obtained, all variables in the hyperbolic system except the radiation energy 
density are updated as usual. For the radiation energy density, two steps are taken. The radiation energy density is 
first updated according to Equation (|45|) with its RHS computed using the Godunov states. Then another Godunov 
solver is used to evolve Equation (|4T|) . which is equivalent to Equation (|4l)|) . The conversion between E g and E v is 
performed according to E g = E u (v g )v g £\\nv g . The Godunov states at time t n+1 / 2 of the first Godunov solver are used 
in computing V u and Vu to obtain the wave speed a q in the frequency space advection. The second Godunov scheme 
for the frequency space advection is a standard approach based on the method of lines. For time integration we use the 



8 



Zhang ct al. 



third-order total variation diminish Runge -Kutta scheme of IShu fe O shcr (1988). We use the piecewise linear method 
with a monotonized central slope limitcr ([van Leenfl977h to reconstruct qi and a q in order to ac hieve high order in 
(logarithmic fre quency) space. The reconstructed variables are used to computed the HLLE flux (|Harten et al.lH983t 
lEisenstatl 1 1 98 lh . The frequency space advection has its own CFL condition for stability. Since CASTRO uses the 
hydrodynamic CFL condition with other additional constraints for the time steps, if necessary, subcycling in time is 
employed in order to satisfy the frequency space CFL condition. 



3.2. Implicit Solver for Radiation Diffusion and Source-Sink Terms 

The implicit solver evolves the radiation and matter according to equations (j4"B"[) . (|4"9")l . and ([50)) with the results of 
the explicit update (S 13. 1| as the initial conditions. Our scheme is based on a first-order backward Euler method. 

3.2.1. Outer Newton Iteration 
We solve the parabolic subsystem (Equations 1481 l49l and [50]) iteratively via Newton's method. We define 

F e = pe - pe~ - At^c^KgEg - jg), (52) 
g 

F Y = P Y e - pY- -AtJ2 c£ ff (%£ 9 - jg), (53) 
g 

F g = E g - Eg- - At V • (^-VE g j + At c{n g E g - j g ), (54) 

where the — superscript denotes the state following the explicit update, and g = 1,2, . . . , N . The desired solution is 
for F g , F e and Fy to all be zero. We approach the solution by Newton iteration, and in each iteration step we solve 
the following system 



" {dF e /dT)W (dF e /dY e )W (dF e /dE g )W 
(dFy/dT)^ (dF Y /dY e )W (dF Y /dE g )W 

_ (OFg/dT)^ (dFg/dY e )W {dFg/dEg)^ 

Here ST^+V = T( fe+1 ' -TW, <5Y"J fe+1) = Y c {k+1) ~Y c {k) and SE g k+1) = E g k+1) - E g k) , where the (k) superscript denotes 
the stage of the Newton iteration. The exact Jacobian matrix in equation ([55[) is very complicated. For simplicity, we 
neglect the derivatives of the diffusion coefficient c\ g /xg- This is a common practice in radiation transport, and does 
not appear to significantly degrade the convergence rate. With the approximate Jacobian matrix, we can eliminate 

ST and SY e from equation (j55[) and obtain the following equation for E g k+1 \ 





- ^(fc+i) " 




- _ F w - 




8Y {k+D 




~ b 9 




_ SE^ _ 





(55) 



CKa + 1_) E (M+i) _ V .f ^V# +1) 



At 



Xg 



c Jg 



El 
At 



+ e 



C E(v4' +1) -v)"Ai^ (fc) -^") 



(56) 



To reduce clutter, we have dropped the (k) superscript for X g , K g , Xgi an d jg- The newly introduced variables in 
equation ([56| arc given by 
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where 



Vy 



By 



cAt 

~~n~ 

cAt 

~~n~ 

cAt 
~~n~ 

cAt 



cAtj2i 



dj__ 

9 \ Q T 



dfi a_ E (k) 



de 



dT 9 
dj dn 



P dYe +C/Xt 2^{ dYe dY, 9 



de 



dT dT ~ 9 



dK 9 E (k) 



n = 



de 



P dT + ^\dT 8T 9 



de 



cAtJ2i 



dj g dn 6 



9\ dT dT a 



(59) 
(60) 
(61) 
(62) 



(63) 



Here all the derivatives are computed from the results of the previous Newton iteration (i.e., the (k) state). 



3.2.2. Inner Iteration 

We now present our algorithm for solving equation (|56|) , which is actually a set of ./V equations coupled through the 
summation terms on the RHS of equation (|5<3|) . The system is usually too large to be solved directly. Instead, we solve 
it by decoupling the groups and utilizing an iterative method. Note that the iteration here for solving equation ()56|) is 
different from the Newton iteration for solving the entire parabolic subsystem. The inner iteration here is embedded 
inside the outer Newton iteration step. In our algorithm, we avoid the coupling by replacing the radiation energy 
density in the coupling terms with its state from the previous inner iteration. Thus, each step of the inner iteration 
solves 



A* 



Xg 



5_ 

At 



e 



g' 

- 3*) - - P Y < 



(64) 



where I is the inner iteration index. Here, we have dropped the (k + 1) superscript for E g +1 and Eg to reduce 
clutter. For the first inner iteration, the state of the previous inner iteration is set to the result of the previous outer 
iteration. The inner iteration for equation is stopped when the maximum of ^2 g \{E g +1 — E g )\/ ^2 g E g +1 on the 

computational domain is below a preset tolerance (e.g., 10~ 6 ). 
CASTRO can solve equations in the following canonical form (Paper II) 

i \ / i i 

where A are cell-centered coefficients and Bi, d, and Di are centered at cell faces. For equation (|64]l the Ci and 
Di coefficients are not used. The same canonical form works for Cartesian, cylindrical, and spherical coordinates 
so long as appropriate metric factors are include d in the coefficients and the RHS. CASTRO uses the hypre library 
(jFalgout fc Yandl2002t \hvvre Code Proiectll2012t ) to solve the linear system. 



3.2.3. Matter Update 
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When the numerical solution to equation (|56l) is obtained via the inner iteration (Equation l64[) . we update the matter 
and continue with the next outer iteration until convergence is achieved. Using Equation (|55[) . we can obtain 



5T (k+i) _ 



E^4 fc+1) - 



■hi ) 



L g 



pe^ — pe 
cAt 



k+1) 



pY™ - pY e 



pe 



k) 



cAt 

- pe~ 



(66) 



cAt 



P Yj k) - P Y- 
cAt 



Thus we could update the matter temperature as T( fe+1 ) = + ST^ 1 ^ and electron fraction as Y"i fe+1 ' = Y, 



(67) 



(fe) 



5Y} k+1 \ and then compute the matter specific internal energy e^ k+1 ^ from T( fc+1 ) and Y t 



(fc+i) 



However, the total 



energy (radiation energy plus matter internal energy) and lcpton number would not be conserved in this type of 
updates. Furthermore, our experience with this approach is that sometimes the update gives T and Y e that are out of 
the bounds of opacity and EOS tables used in our core-collapse simulations and hence causes convergence issues. Thus 
a different approach is taken in CASTRO. To conserve the total energy and lepton number, we update the matter as 
follows, 



(fe+i) _ 



pe 



?Yj k+1 



Hpe [ 
+ cAt 



f Q{pY^-pY-) 
j g )-(H + QZ g ){ Kg El 



+ (1 - H)pe 

E l^ E t +1 

a 

QpY^ + (1 - Q)pY e - + H(peW - pe ~) 
+ cAtY J [U* g El +1 



■j a )} 



(68) 



(69) 



where H = ^2 g H g , O = J2 g ®g > H = S g €gHg> an d © = ^2 g £g®g- Then we compute temperature T^ k+1 ^ and electron 

fraction Ye k from the newly updated pe ( - fc+1 - ) and pYe k+1 \ We also compute the new absorption, scattering and 
emission coefficients using the new temperature and electron fraction. By updating the matter using equations (|68[) 
and (I69|) . the implicit algorithm conserves the total energy and electron lepton number regardless of the accuracy of 
the solution to the parabolic subsystem. 
The outer iteration for solving equations (|48[) - (|50l) uses the following convergence criteria: 



\ T (k+i) 

|y e (fc+l) 



_ j.(fe)i 

-y e ( fc )| 



< eT (k+i) 

<e y e (fc+l) 

de 



\F^ 



<^|l(^) (fc+1) ll^ (fc+1) -r (fe) | + |(^) ( ' 
<epYj k+1 \ 



(70) 
(71) 

(72) 
(73) 



where e is a small number such as e = 10~ 6 . We do not use the relative error of the matter specific internal energy e 
because it may not accurately reflect the error since that e can be shifted by an arbitrary amount. Typically, the outer 
iteration is stopped only when all four conditions ( Inequalities 1701 - I73|) are satisfied. But for problems with extremely 
high opacities, Inequalities (|72[) & (|73]l can be difficult to achieve, and only Inequalities (|70[) & (|7"Tj) are used. 

3.2.4. Acceleration Scheme 

We solve the multigroup diffusion equation (Equation I56[) using an iterative method. The system of the multigroup 
radiation diffusion equations is coupled. Each individual radiation group depends on the matter, which in turn depends 
on all radiation groups. When the coupling is strong (e.g., large opacity and s mall he at capacity) and/or the time 
step is large, the convergence rate of the iterative method we have presented in jjl3~2.2l ca n be v ery slow. To improve 
the convergence rate, we have extended the synthetic acceleration scheme of iMorel et al.l (|1985| ) for photon radiation 
to neutrino radiation, and have developed two types of acceleration schemes. 



Local Acceleration Scheme — We define the error of radiation energy density after £ inner iterations as 

d +1 = Et - Ei +1 , 



(74) 
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where E g is the exact solution of Equation (|56p and E g +1 is the result of the i-th inner iteration via Equation (|M|) . 
In our local acceleration scheme, the spatial derivatives of e e g +1 are assumed to be zero. This is justified because the 
spatially constant limit corresponds to the most slowly convergent mode ([Morel et al.lll985l ) and our goal is to improve 
the convergence rate by curing the slowest decaying mode of the errors. Using equations (|56[) and (|64p . we derive an 
equation for the error in the g-th group, 

44 +1 - H a E - 6 » £ = Brf 1 + ^'y\ (75) 

a' a' 

where 



«* a = "9+z^l, (76) 



r t+1 

T 



r e+1 
'y 



Ev(^-4). ( 77 ) 



E^(^ +1 -4)- (78) 

9' 



Eq. [75] can be solved analytically, and gives 

e l+1 = 

9 Kgipq-rs) 



e i+i = (gg g + r9 g )4 +1 + ( P Q g + sH g y Y +1 ^ (7g) 



where 



?=i-E% (81) 



9 

£^r^ (82) 



9 



The improved solution is therefore 



E%?- (83) 



£^+ 1 :=^+ 1 + 4+ 1 . (84) 



The acceleration is placed after the convergence check for the inner iteration so that no unnecessary acceleration is 
applied to converged solutions. 

Gray Acceleration Scheme — The original acceleration scheme of lMorel et al.l (jl985|) utilizes a gray diffusion equation for 
errors. It is postulated that the normalized multigroup spectrum of the errors is given by the normalized equilibrium 
spectrum. Note that our local acceleration scheme is based on the assumption of local equilibrium. Therefore the 
spectrum defined as 



^ g = ;^4+T- (85) 

^9' V 



is postulated to be 



c 



■e+i 

9 



^9 = ^ Vl ' ( 86 ) 



a' 9' 



where e g +1 is the solution of the local acceleration scheme (Equation [79]) . Substituting E g +1 = E g — e g +1 into Eq. [64] 
and summing over groups, we obtain 

Ae e+1 - V ■ (dVe £+1 ) + V ■ (0e £+1 ) = cHr e T +1 + cOr e Y +1 , (87) 
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where 

a 

A = c(l — H)R - c<dR Y + (89) 

d ~ = E^> ( 9 °) 

^=-E^ V ^' (91) 

s = E^9 K S' (92) 
s 

Ky = ^^gigKg. (93) 

Equation (|87|) is in the canonical form that CASTRO solves (Equation [53]). The solution of the gray diffusion equation 
for error, e +1 , is used to obtain the improved multigroup solution 

E e g +1 :=^+ 1 +^+ 1 . (94) 

We note again that the acceleration operation is placed after the convergence check for the inner iteration so that no 
unnecessary acceleration is applied to converged solutions. 

4. TEST PROBLEMS 

In this section we present detailed tests of the code demonstrating its ability to handle a wide range of radiation 
hydrodynamics problems. Photon radiation problems are presented in Section [4. 1\ whereas Section f4. 21 shows neutrino 
radiation hydrodynamics simulations of core-collapse supernovae. Note that not every term or equation in the full 
system of radiation hydrodynamics is included in every test. This allows various parts of the code to be tested 
separately. 

A CFL number of 0.8 is used for these tests unless stated otherwise or a fixed time step is used. For photon radiation 
test problems, no additional constraint is posed by the implicit solver. For core-collapse supernova simulations, the 
time step is further constrained by limiting the estimated relative changes in temperature to less than 1% and the 
estimated absolute changes in electron fraction to less than 0.01 using the change rates from the previous time step, and 
in our simulations this only happens around the time of core bounce. The relative tolerances for both the outer and the 
inner iterations in the implicit update are 10~ 6 . All convergence criteria for the outer iteration (Inequalities 1701 - IT3"f 
are used except for the tests in Sections 14.1 . 41 and 14 . 1 . 71 with extremely high opacities. By default, the radiation groups 
are uniform in logarithmic frequency space, and the representative frequency for the g-th group is v g = ^/v g -i/2 v g+i/2- 
For AMR simulations, the "deferred sync" scheme (Paper II) is used for flux synchronization between coarse and fine 
levels in the implicit solver. 

4.1. Photon Radiation Test Problems 
4.1.1. Linear Multigroup Diffusion 

In this problem, we simulate a linear multigroup diffusion test problem introduced by iShestakov fc Bolstadl (|2005f) . 
This can test the multigroup diffusion part of the system 

3 

= " <^ E 9 3„) + V • (^ E a) ■ (96) 

To make the nonlinear system mathematically tractable, [Shcstakov & Bolstad ( 2005|) linearize the system using the 
following assumptions. First, the absorption coefficient is assumed to be proportional to and the group-averaged 
absorption is given by 

k« = C K v-\ (97) 



where C K is a constant. Second, radiation emission is assumed to obey Wien's law rather than Planck's law, and it is 
further modified so that 

hv a -i/2\ ( hv g+1 / 2 " 



B a = —^a 



cxp — — — cxp 



kT f r V kT< 



T, (98) 
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Fig. 1. — Linear multigroup diffusion test. Numerical results are shown in lines, whereas the analytic solutions are shown in circle 
symbols. We show dimcnsionlcss total radiation energy density and temperature at the dimensionless time of 1 in panels (a) and (b), 
respectively. The spatial coordinate x is also dimensionless. 



TABLE 1 

Relative errors of linear multigroup diffusion test 
problem. The relative errors of temperature and total 
radiation energy density at the dimensionless time of 1 are 

SHOWN. 
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e(E r ) X 10 3 


0.00 


0.0017 
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4.8283 
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0.20 
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0.40 
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0.46 
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0.54 
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0.4443 
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0.0461 
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0.49 
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1.00 


1.3612 


2.2413 


0.50 


0.0020 


0.3712 









a Dimensionless spatial coordinate 



where Tf is a fixed temperature. Note that the group-integrated B g now has linea r dependence on T. F i nally, there 
is no scattering, and the diffusion is not flux limited. Under these assumptions, iShestakov fc Bolsta d (2005) have 

obtained exact solutions for the linearized system. 

The test we have simulated is the first prob lem in[S hestakov & Bolstad (2005]). We have set up the test using physical 
units rather than the dimensionless units of IShestakov &: Offnerl (|2008|) . The one-dimensional computational domain 
in physical units is (0,1.0285926 x 10 6 ) cm, corresponding to (0,5.12) in dimensionless units. We use a symmetry 
boundary condition at the lower boundary, and a Marshak boundary with no incoming flux at the upper boundary 
(i.e., Eg + (2/3n g )dE g /dx = 0). The mass density is 1.8212111 x 10~ 5 gcm~ 3 . The specific heat capacity of the 
matter at constant volume is 9.9968637 x 10 7 ergK^ 1 g _1 . Initially, the temperature is set to 0.1 keV and 0, for 
x < 0.5 and x > 0.5, respectively. Here, a; is a dimensionless coordinate. The radiation energy density is initially 
set to zero everywhere. The constant C K is set to 4.0628337 x 10 43 cm _1 Hz 3 . The fixed temperature in the group- 
integrated radiation emission (Equation l98|) is Tf = 0.01 keV, corresponding to Tf = 0.1 in dimensionless units. We 
use 64 radiation groups with the lowest boundary at zero. The width of the first group is set to 1.2089946 x 10 13 Hz, 
corresponding to 5 x 10 -4 in dimensionless units, and the width of other groups is set to be 1.1 times the width of its 
immediately preceding group. The representative frequency for each group is set to u g = ^/^ 3 -i/2 fg+1/21 except that 
v\ is set to 0.5^3/2 for the first group. We have run a simulation with 2048 uniform cells. Thus the dimensionless 
cell size is 1/400. The time step is fixed at 5.8034112 x 10 _8 s, corresponding to 1/200 in dimensionless units. No 
acceler ation scheme is used in this test. Figure Q] shows the results after 200 steps. The comparison with the tabular 
data of IShestakov fc Off ncr (2008|) is shown in Tablc[T] Relative errors arc computed as e(/) = |/ n — / |// c , where / n 
and f e are the numerical and exact results, respectively. Because the tabular data of the exact solution are located 
on the faces of numerical cells, averaging of the numerical data has been performed for the com parison. The results 
are in excellent agreement with the exact solution, and are comparable to the numerical results of IShestakov fc Offnerl 
(|200ll . 
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Fig. 2. — Noncquilibrium radiative transfer with picket-fence model. Numerical results of Ui, U2, and V are shown for t = 3 (solid 
lines) and 30 (dashed lines). Also shown are analytic results for r = 3 (circles) and 30 (diamonds). 

4.1.2. Nonequilibrium Radiative Transfer With Picket-Fence Model 

ISu fc Olsonl (|1999f ) have studied a nonequilibrium radiative transfer problem with a multigroup model, more specif- 
ically the picket-fence model. There are some special assumptions in this test problem. The volumetric heat capacity 
at constant volume is assumed to be c v = aT 3 , where T is the matter temperature and a is a parameter. The 
group-integrated Planck function is assumed to be 



(99) 



where a is the radiation constant, and p g are parameters such that ^2 Q P g = 1- It is also assumed that the absorption 
coefficient is independent of temperature and there is no scattering. ISu fc Olsonl (|1999D have also defined dimensionless 
variables for convenience. The dimensionless spatial coordinate is defined as 

(100) 



KZ, 



where z is the coordinate in physical units, and 



The dimensionless time is defined as 



g 



t. 



(101) 



(102) 



where t is the time in physical units. The dimensionless variables U g for radiation energy density and V for matter 
energy density are defined as 



V 



(103) 
(104) 



where To is a reference temperature. In this test problem, t here is a con s tant r adiation source that exists for a finite 
period of time (0 < r < To) on a finite range (|x| < xo). ISu fc Olsonl ()1999l ) have found analytic solutions to the 
problem for both transport and diff usion approaches. 

Here, we simulate the Case C of ISu fc Olsonl (|1999( ) by solvi ng the mult i group diffusion equations (Equations 1951 
& in!]), and compare the results with the diffusion solution of ISu fc Olsonl (p~999f ) . Two radiation groups are used 
in this test, with pi = P2 = 0.5. The absorption coefficients are set to k\ = 2/101 cm -1 and K2 = 200/101 cm" 1 . 
Thus, R = lcm -1 . The constant a in the expression for volumetric heat capacity is set to a = 4a. The reference 
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Fig. 3. — Spectrum from the radiating sphere observed at t = 10 12 cm and r = 0.06 cm. Both numerical (circle symbols) and analytic 
(solid line) results are shown. 

temperature is chosen as To = 10 6 K. The parameters for the radiation source are xq = 0.5, and tq = 10. When r < To, 
the radiation source deposits radiation energy to the region |x| < xq and increases the radiation energy density at a 
rate of PgCRaT^ for each group. The computational domain is < x < 102.4 covered by 1024 uniform cells. The lower 
boundary is symmetric, and the radiation energy density outside the upper boundary is set to zero. Initially, there is 
no radiation and the temperature is zero. A fixed time step At — 0.1 is employed. The explicit solver is turned off 
in the simulation, and the radiation flux is not limited. No acceleration scheme is used in this test. The numerical 
results are shown in Figure [2j Our results are in excellent agreement with the analytic solution. 

4.1.3. Radiating Sphere 

Graziani (2008) found an analytic solut ion for the spec t rum o f a radiating sphere in a medium with frequency- 
dependent opacity. This test was used by iSwestv fc Mvral (|2009f ) to test their V2D code. In this test, a hot sphere 
with a fixed temperature is surrounded by a cold static medium. The absorption coefficient is assumed to be zero, 
and the radiative transfer is due to coherent and isotropic scattering. Thus, there is no radiation-matter coupling or 
group coupling. The time-dependent spectrum at a given place in the medium is therefore the superposition of the 
black-body spectrum of the medium and the spectrum of the radiation originated from the hot sphere. 

We perform this test in one-dimensional spherical geometry. The computational domain is 0.02 cm < r < 0.1 cm and 
is covered by 128 uniform cells. The hot sphere has a fixed temperature of 1.5 kcV, and the medium is at 50 eV. We 
use 60 energy groups covering an energy range of 0.5 eV to 308 keV. The group Rosseland mean coefficient is assumed 
to be Xg = 10 13 (z/ s /3.6 x 10 14 Hz) -3 cm -1 . The explicit solver is turned off. No acceleration scheme is used in this 
test. There is no outer iteration in the implicit solver because the matter properties do not change. The time step is 
fixed at 10 _15 s. The spectrum at r — 0.06 cm and t = 10~ 12 s is shown in Figure |3] The numerical results are again 
in good agreement with the analytic solution. 

4.1.4. Shock Tube Problem In Strong Equilibrium Regime 

In this one-dimensional problem, the initial conditions of the matter arc p^ = 10~ 5 gcm~ 3 , Tl = 1.5 x 10 6 K, 
ul = and pr = 10 _5 gcm -3 , Tr = 3 x 10 5 K, ur — 0, where L stands for the left half, and R the right half of 
the computational domain (0 < x < 100 cm). The gas is assumed to be ideal with an adiabatic index of 7 = 4/3 
and a mean molecular weight of /1 = 1. Initially, the radiation is assumed to be in thermal equilibrium with the gas 
(i.e., Eg — AnBg/c). The interaction coefficients arc set to K g = 10 6 cm _1 and Xg = 10 s cm -1 . Thus, due to the 
huge opacities, the system is close to the limit of strong equilibrium with no diffusion, and is essentially governed 
by Equations (|3"U1) (|3~4"|) with A g = 1/3 and f g = 1/3. The parameters of this test problem are chosen such that 
neither radiation pressure nor gas pressure can be ignored. This test is adapted from the shock tube problem in Paper 
II. However, the adiabatic index of the gas is now 7 = 4/3 rather than 5/3 so that we can compare the numerical 
results with the exact solution of a pure hydrodynamic Riemann problem. Furthermore, although this is a one- 
temperature gray problem due to the huge opacity, the numerical calculation is performed with 16 radiation groups 
with the lowest and uppermost boundaries at 10 14 Hz and 10 18 Hz, respectively. We use 128 uniform zones in this test. 
The numerical results are shown in Figure 2] Also shown are the exact solutions of a pure hydrodynamic Riemann 
problem corresponding to the numerical problem. The results show good agreement. We also not e that our scheme is 
stable without using small time steps even though the system is in the "dynamic diffusion" limit (lMihalas &: Mihalai 
Il999t) with tu/c ~ 10 7 , where r is the optical depth of the system. Either the gray acceleration scheme or the local 
acceleration scheme (§ 13.2. 4[) must be employed in this test, otherwise the implicit solver, which is responsible for the 
thermal equilibrium between radiation and matter, exhibits extremely slow convergence. 
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Fig. 4. — Shock tube at t = 10~ 6 s. Numerical results from a full radiation hydrodynamics calculation with 16 radiation groups and 
128 cells are shown in symbols. The exact solutions of a pure hydrodynamic Ricmann problem corresponding to the numerical problem 
are shown in solid lines for comparison. We show mass density (p), velocity (v), total pressure (ptot), and radiation energy density (E r ). 
For the radiation hydrodynamics simulation, the total pressure is the sum of gas pressure and radiation pressure, and the radiation energy 
density is summed over all groups. In the pure hydrodynamic Ricmann problem, there is no radiation energy density. All quantities are in 
cgs units. 
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Fig. 5. — Density and temperature profiles for Mach 5 supercritical shock. Numerical results are shown in symbols. We show density 
{circles), gas temperature (plus signs), and radiation temperature (squares). Here radiation temperature is defined as (Er/a) 1 ^, where E r 
is the total radiation energy density. Also shown are the analytic results for density (red solid line), gas temperature (blue solid line), and 
radiation temperature (green solid line). We show only part of the region near the density jump. The inset shows a blow-up of the spike 
in temperature. The numerical results have been shifted by —207 cm in space to compensate for the discrepancy in shock position caused 
by the initial setup. 

4.1.5. Nonequilibrium Radiative Shock 

In this test, we simulate the Mach 5 shock problem of lLowrie k, E dwards (20Q8[). In Paper II, we showed results 
computed with the gray radiation hydrodynamics solver. Here, we present the numerical results obtained with the 
multigroup radiation hydrodynamics solver. In this test, the one-dimensional numerical region (—4000 cm, 2000 cm) 
initially consists of two constant states: pr, = 5.45969 x 10 -13 gcm -3 , Tl — 100 K, mj, = 5.88588 x 10 5 cms _1 and 
Pr = 1.96435 x 10- 12 gcuT 3 , T R = 855.720 K, u R = 1.63592 x lO^ms" 1 , where L stands for the left, and R the 
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right of the point x = 0, respectively. The gas is assumed to be ideal with an adiabatic index of 7 = 5/3 and a mean 
molecular weight of \i = 1. The setup is the same as that in Paper 110 except that 16 radiation groups are used in 
this paper. The lowest and uppermost boundaries of the groups are at 10 10 Hz, and 10 15 Hz, respectively. Initially, the 
radiation is assumed to be in thermal equilibrium with the gas (i.e., E g = 4irB g /c). The interaction coefficients are 
set to K g = 3.92664 x 10~ 5 cm -1 and Xg = 0.848903 cm -1 , respectively. The states at spatial boundaries are held at 
fixed values. Four refinement levels (five total levels) are used with a refinement factor of 2 with the finest resolution 
being Ax rj 1.5 cm. In this test, a cell is tagged for refinement if the second derivative of either density or temperature 
exceeds certain thresholds (Paper II). The initial conditions are set according to the pre-shock and post-shock states 
of the Mach 5 shock. After a brief period of adjustment, the system settles down to a steady shock. The simulation 
is stopped at t = 0.04 s. The results are shown in Figure [SJ The agreement with the analytic solution including the 
narrow spike in temperature is excellent, except that the numerical results in the figure had to be shifted by —207 cm 
in space to match the analytic shock position. This discrepancy in shock position is due to the initial transient phase 
as the initial state relaxes to the corre c t stea dy-state profile. No flux limiter is used in this calculation because the 
analytic solution of iLowrie fc Edwards! (pOOSh assumes A = 1/3. We again found that the implicit solver would fail 
to converge without an acceleration scheme. Both the gray acceleration scheme and the local acceleration scheme 
(§ ET2T4]) work in this test. 

4.1.6. Advecting Radiation Pulse 

In Paper II, we performed a test of advecting a ra diation pulse by the motion of matter with the gray radiation 
hydrodynamics solver (see also iKrumholz et al.l 12007). Here, we repeat the test with small changes to exercise the 
multigroup aspect of CASTRO. The setup is the same as that in Paper II, except that the absorption coefficient k in 
this paper is given by 
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As in the test in Paper II, the initial temperature and density profiles are 



cm" 1 . (105) 



r = r + (T 1 -To)cxp(-— ), (106) 
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"»5T + S(f- TS j- < 107 » 

where To = 10 7 K, T\ = 2 x 10 7 K, po = 1.2 gem -3 , w = 24cm, the mean molecular weight of the gas is /x = 2.33, 
and R is the ideal gas constant. Initially the radiation is assumed to be in thermal equilibrium with the gas. No flux 
limiter is used in this test (i.e., A = 1/3 and / = 1/3). If there were no radiation diffusion, the system would be in 
an equilibrium with the gas pressure and radiation pressure balancing each other. Because of radiation diffusion, the 
balance is lost and the gas moves. 

The purpose of this test is twofold. First, we use this test to assess the ability of the code to handle radiation 
being advected by the motion of matter. We solve the problem numerically in different laboratory frames. Then 
we compare the case in which the gas is initially at rest to the case in which the gas initially moves at a constant 
velocity. Another purpose is to assess the effect of group resolution. Since the opacity is very large, the system is 
close to thermal equilibrium between radiation and matter. Hence, the results of multigroup radiation hydrodynamics 
simulations should become increasingly close to those of a gray radiation hydrodynamics simulation as more groups 
are used. 

We have performed four runs: an unadvected run with 8 radiation groups, an advected run with 8 groups, an advected 
run with 64 groups, and an advected run using the gray radiation hydrodynamics solver presented in Paper II. The 
velocity in the unadvected run is initially zero everywhere, whereas in the advected runs it is 10 6 cms _1 everywhere. 
For the multigroup runs, the lowest and uppermost group boundaries are at 10 15 Hz and 10 19 Hz, respectively. In the 
gray simulation, the single group Planck mean coefficient, Kp, and Rosseland mean coefficient, kr , of the opacity given 
by Equation (|105j) are 

k p = 3063.96 (-^—] cm" 1 , (108) 




k r = 101.248 (j^) cm- 1 . (109) 

The computational domain in all four runs is a one-dimensional region of —512 cm < x < 512 cm with periodic 
boundaries. The numerical grid consists of 512 uniform cells. The simulations are stopped at t = 4.8 x 10 ~ 5 s. 
Figure |JJ] shows the density, velocity, and temperature profiles for all four runs. The results of the advected runs 
have been shifted in space for comparison. The profiles from the unadvected 8-group run and the advected 8-group 



4 Unlike the simulation shown in Paper II, the simulation here used the 2010 CODATA recommended values of the fundamental physical 
constants. Therefore, the physical values of the states have changed slightly. 



18 



Zhang et al. 



1.2 
^ 1.1 

\ 10 

°- 0.8 
0.7 

6 
4 

T » 2 

a o 

==- -2 
3 -4 
-6 

1.30 
1.25 
* 120 

s 1,15 
g-i.io 

1.05 
1.00 




— unadvected, 8 grps 

- - advected, 8 grps 



— advected, gray 

- - advected, 64 grps 




Fig. 6. — Density, velocity and temperature profiles at t = 4.8 X 10 s for the test of advecting radiation pulse. We show the results 
of four runs including the unadvected 8-group run (red solid lines), the advected 8-group run (blue dotted lines), the advected 64-group 
run (black dotted lines), and the advected gray run (green solid lines). The results of the advected runs have been shifted in space for 
comparison. 
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Fig. 7. — Relative differences in density (solid lines) and gas temperature (dashed lines) in the test of advecting radiation pulse. The 
relative differences between the advected 8-group and unadvected 8-group runs are shown in panel (a), whereas the relative differences 
between the advected 64-group and advected gray runs in panel (b). 

run are almost identical demonstrating the accuracy of our scheme in handling radiation advection by the motion of 
matter. This is expected because the velocity is significantly smaller than the speed of light (v/c ~ 3 x 10~ 5 ). The 
results from the advected 64-group run and the gray run are also almost identical as expected. It is clear that the 
8-group results differ from the gray results. This is not surprising because we did not perform group averaging for 
interaction coefficients (Equations [21] & [22]) . We also show the relative differences in Figure \7\ The relative difference 
of the advected 8-group run with respect to the unadvected 8-group run is computed as (advected - unadvected) / 
unadvected, and the relative difference of the 64-group run with respect to the gray run is computed as (multigroup - 
gray) / gray. The figure shows that the difference between the advected 8-group and unadvected 8-group runs is less 
than 0.03% everywhere, and the difference between the 64-group and gray runs is less than 0.1% everywhere. 

4.1.7. Static Equilibrium 

In Paper II, we have demonstrated that the gray solver can maintain a perfect static equilibrium in multiple dimen- 
sions because the radiation pressure and the gas pressure are coupled in the Riemann solver. Here, we repeat the test 
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Fig. 8. — Structure of the magnitude of velocity at t = 10 -4 s in the static equilibrium test. The unit of velocity is cm s _1 . 

using the multigroup solver. The setup is the same as that in Paper II. The initial conditions for radiation and matter 
are the same as those in the test of the advecting radiation pulse (§ 14. 1 .6[) except that the interaction coefficients are 
set to a very high value, n g = \ g = 10 cm g x p. Thus, almost no diffusion can happen, and the radiation and 
the matter are in perfect thermal equilibrium. We have performed a calculation on a two-dimensional Cartesian grid 
of —512 cm < x < 512 cm and —512 cm < y < 512 cm with 512 uniform cells in each direction. The initial velocity is 
zero everywhere. For the initial setup, the coordinate x in equation (|106[) is replaced by r = \J x 2 + y 2 . We use 16 
radiation groups with the lowest and uppermost boundaries at 10 15 Hz and 10 19 Hz, respectively. Figure [8] shows the 
velocity profile at t = 10~ 4 s (after about 6000 steps). The maximal velocity at that time is ~ 8 x 10~ 8 cm s _1 , which 
is about 10~ 15 of the sound speed. Such a small gas velocity indicates that the multigroup solver in CASTRO can 
also maintain a perfect static equilibrium in multiple dimensions, as does the gray solver in CASTRO. 

4.2. Core- Collapse Supernova Simulations 

In this section, we present core-collapse supernova simulations in ID spherical and 2D cylindrical coordinates. 
Newtonian gravity with a monopole approximation is used in the 2D simulation for simplicity, although CASTRO 
has a multigrid Poisson solver for gravity. In the simulations, there are three neutrino species: electron neutrinos, 
v e , electron antincutrinos, v e , and which stands for a combined species representing i/p, v T , Dp, and v r . We 
use a tabular equation of state that provides thermodynamic qu antities as a function of temperature, density, and 
electron fraction and was constructed from the iShen et al.l (1998a. b) mean-field table, augmented with photons and a 
ge neral electron /positro n equation of state, and extended down to 10 g cm~ 3 . The opacities used arc fully described 
in lBurrows et al.l (|2006Q . They include 1) elastic scattering off nucleon s, alpha par t icles, a nd the single representative 
heavy nucleus at the peak of the binding energy curve calculated by IShen et all (11998alTQ), a nd 2) charged-current 
absorption cross sections off nucleons and nuclei (the latter using the approach of lBruennl (fl985D ). The scattering and 
absorption cross sections off nucleons are corrected for weak magnetism and recoil effects, the form factor, electron 
screening, and ion-ion correlation effects on Freedman scattering off nuclei are applied, and pair production by nucleon- 
nucleon bremsstrahlung and electron-positron annihilation (and their inverses) are included as sources (and sinks). 
The sink terms of the latter do not include Fermi blocking corrections by the neutrinos in the final state. Inelastic 
neutrino-electron scattering is treated as elastic scattering. 

4.2.1. One- Dimensional Core-Collapse Supernova Simulation 

We present here a core-colla pse supernova simulation in one-dimensional spherical coordinates. The initial model is 
based on a 15 M Q progenitor of lWooslev k, Weaver! (|1995[ ). The computational domain for this run is < r < 5x 10 8 cm. 
The simulation employs four refinement levels (five total levels) with 1024 zones on the coarsest level (i.e., level 0) and 
a resolution of Ar ps 0.305 km on the finest level (i.e., level 4). The AMR refinement factor is 2. The regions where the 
density is greater than 7.2 x 10 6 g cm -3 , or the enclosed mass is greater than 1.32 M Q are refined to level 2. In addition, 
cells on level I will be refined to level £+1 if Y e < 0.4, s > 5, VF e > 0.01/Ar^, or Vs > 1/Are, where I < 4, s is entropy 
in units of /cbaryon -1 , and Are is the cell width on level £. We use 40 energy groups for electron neutrinos with the 
first group at 1 McV and the last group at 300 MeV. We use 32 energy groups for electron antineutrinos with the first 
and last groups at 1 MeV and 100 MeV, respectively. For v a , we use 40 energy groups with the first and last groups at 
1 McV and 300 McV, respectively. The flux limiter of iLevermore &: Pomrani ng (1981j) is used is this simulation. The 
local acceleration scheme is employed. The inner boundary is symmetric. An outflow boundary condition is used for 
the outer boundary in the explicit solver. The modified Marshak boundary of iSanchez fc Pomraning] (|1991|) with zero 
incoming flux is used for the outer boundary in the implicit solver. Thus, the region outside the outer boundary is 
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Fig. 9. — Density (upper left), velocity (lower left), temperature (upper center), entropy per baryon (lower center), electron fraction 
(upper right), and net lepton fraction (lower right) as a function of enclosed mass for the ID core-collapse supernova simulation. Snapshots 
are shown for (black solid lines), 1 ms (red dashed lines), 15 ms (green dotted lines), and 115 ms (blue dash-dotted lines) after bounce. 

essentially treated as a vacuum, and the radiation diffuses into the vacuum. Note that in the usual Marshak boundary 
based on the classical diffusion theory, the flux at the boundary does not have the correct behavior in the streaming 
limit, and the issue has been fixed in the modified Marshak boundary of iSanchez fe Pomrani ng (1991) by using flux 
limited diffusion. 

The simulation is carried out to 200 ms after core bounce. Figure [9] shows four snapshots of density, velocity, 
temperature, entropy, electron fraction, and net electron lepton fraction (Yt). The infall of the material due to 
gravity results in the increase of the central density until the core bounces. At bounce, the density reaches p corG sa 
2.9 x 10 14 g cm -3 at the core. A shock appears around bounce at M sn « 0.635M Q . The electron and lepton fractions 
at bounce are Y e corc py 0.288 and Y£ oro pa 0.326, respectively, at the core. The core temperature at bounce is T coro sa 
9.98 MeV. T he core entropy at bounce is a core sa 0.901 fc baryon" 1 . These propertie s are in agreement w ith previous 
studie s (e.g., iRampp fc J ankal 120001 : iLiebendorfer et aH 1200 it iThompson et al.ll2003| ). In a recent study, iLentz et"aTI 
(12011) have obtained M sh w 0.492 M Q , p cole w 4.264 x 10 14 gcm" 3 , F e core Pa 0.2407, and Y£ mo Pa 0.2782 for their model 
N-FullOp at bounce, and M sh sa 0.717 M Q , p core sa 3.336 x 10 14 gem" 3 , r e corc sa 0.3046, and Y[ OIC Pa 0.3696 for their 
model N-RcduceOp at bounce. Except for p core , the properties of our model at bounce are in between those in models 
N-FullOp and N-ReduceOp of ILentz et all (|2012j) . 

Figure [9] shows a spike near the shock in the Y e and Yl profiles shortly after bounce (see also ILentz et al.ll2012ft . 
The spike of Yl is caused by electron neutrinos leaking out of the shock. The absorption of some of those neutrinos 
by the material in front of the shock results in a spike of Y e . Behind the shock, the minimal Y e drops rapidly for the 
first few milliseconds after bounce due to electron capture; it then becomes steady at ~ 0.065 — 0.075. The velocity 
in some region behind the shock is positive for the first ~ 2 ms after bounce; th e max imum positive velocity reaches 



2 x 10 9 cm s _1 similar to the maximum positive velocity in IThompson et al.l ([2003D . The maximum temperature 
in the shocked region increases rapidly to 17.5 MeV at ~ 0.1 ms after bounce; it then decreases quickly as the shock 
expands outward in radius; it increases again and reaches ~ 22 MeV at 200 ms after bounce due to compression and 
deleptonization as matter settles onto the proto-neutron star. The evolution of the temper ature, electron fraction and 
velocity profiles is in agreement with that in previous studies (e.g.. IThompson et al.l[2003l ). 

In Figure I1Q[ we show shock radius as a function of time. For the first few milliseconds after bounce, the shock 
expands very rapidly in radius, and reaches ~ 72 km. The peak shock radius reaches ~ 140 km at ~ 100 ms after 
bounce. After that the shock slowly shrinks in radius although it continues to grow in mass. 

The comoving frame luminosities of v e , v e , and v u neutrinos measured at 500 km are shown i n Figure 111! In 
agreement with previous studies (e.g.. IThompson et alJl2003t iMarek fc Jankall2009t ILentz et al.ll2012t ). there is a small 
dip in L Ve after the initial rise. The dip is followed by a large pulse in L Vc that reaches 540 Bethe s _1 at 4 ms after 
bounce. The full width at ha lf maximum for L v is ~ 7 ms. The peak L Vc in our model is about 20% larger than 



that in IMarek fc Jankal (|2009l); ILentz et al.1 (|2012j ). In our results, there is also a noticeable spike in L v . Note that 
IThompson et al.l (|2003t ); ILentz et al.l (|2012f) have obtained somewhat similar spikes in their simulations that did not 
include inelastic scattering. The lack of inelastic scattering of neutrinos on electrons results in harder neutrino spectra, 
especially for v^, as shown in the rms average energies (Figure IT2"j) . Here the rms average is defined as 



J n u Av 



1/2 



(110) 



where e v is the neutrino energy and n„ is the monochromatic neutrino number density. Nevertheless we obtain 
(e„ e ) < (e Pc ) < {e v ) after the initial phases, as expected. At 200 ms after bounce, we obtain (e Ve ) sa 14, (e^J sa 19, 
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Fig. 10. — Shock radius vs. post-bounce time for the ID core-collapse supernova simulation. The shock position is defined as the place 
with the largest negative gradient of radial velocity. The peak shock radius reaches ~ 140 km at ~ 100 ms after bounce. 
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Fig. 11. — Comoving frame luminosities of v e (red solid line), u e (green dotted line), and neutrinos (blue dashed line) measured at 
500 km for the ID core-collapse supernova simulation. The luminosity of electron neutrinos reaches a peak value of 540 Bcthc s~ x at 4 ms 
after bounce. For L Vll , the full width at half maximum is ~ 7 ms. There is also a noticeable spike in the luminosity of neutrinos with a 
peak value of 160 Bethe s — 1 at 17 ms after bounce. 
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Fig. 12. — Rms average energies of i/ e (red solid line), v e (green dotted line), and Up neutrinos (blue dashed line) measured at 500 km for 
the ID core-collapse supernova simulation. For (t v ), the peak value of 61MeV is reached at 13 ms after bounce. The rms average energies 
at 200 ms after bounce are 14, 19, and 26MeV, for u e , u e , and v^, respectively. We do not show (e^) before bounce because there is very 
little Ufj, at those early times. 
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Fig. 13. — Luminosity spectrum of u e (red solid line), P e (green dotted line), and neutrinos (blue dashed line) measured at 500 km 
and 200 ms after bounce for the ID core-collapse supernova simulation. The actual energy groups are shown in circles. 
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Fig. 14. — Snapshot at 1.5 ms after bounce for the 2D cylindrical core-collapse supernova simulation. We show (a) logarithmic density, 
(b) entropy per baryon, (c) temperature, and (d) electron fraction, where the unit for density is g cm~ 3 . Only the inner 100 km of the 
model is shown here. 

and (e„ ti ) w 26 MeV. The luminosity spectra measured at 500 km in the comoving are shown Figure [T2J The spectra 
are evidently harder than those in iThompson et al.l (|2003l ) because inelastic scattering is not included in our model. 

4.2.2. Two-Dimensional Core-Collapse Supernova Simulation 

Here we present a two-dimensional core-collapse supernova simulation. The computational domain for this 2D 
cylindrical run is < r < 2.5 x 10 s cm and —2.5 x 10 s cm < z < 2.5 x 10 8 cm. Two refinement levels (three total 
levels) with a refinement factor of 4 are used. The coarsest level has 256 and 512 cells at r-direction and ^-direction, 
respectively. Thus the size of the cells on the finest level is about 0.61 km. The same refinement criteria as those in 
the ID simulation in Section 14.2.11 arc used. We use 20, 16, and 20 groups with logarithmic spacing for z/ e , v e: and v^, 
respectively. The other aspects of this 2D run such as the outer boundary and flux limiter are the same as those in the 
ID simulation in Scction l4.2.1l We start the 2D simulation from the result of a ID simulation rather than the 15 M Q 
prcsupcrnova model. The initial model is obtained from the ID simulation at about 10 ms before bounce. At that 
time, multi-dimensional effects are expected to still be small. The ID simulation that provides the initial model for 
the 2D run has similar resolution as the 2D run. Note that, other than resolution, the low-resolution ID simulation 
here is essentially the same as the ID high-resolution simulation in Section 14.2.11 

Figure fT4l shows the density, entropy, temperature, and electron fraction profiles at 1.5 ms after bounce. The profiles 
are still highly symmetric. Spherically averaged profiles of density, velocity, temperature, entropy, electron fraction, 
and lepton fraction at 1.5 ms are shown in Figure [T5l Also shown in Figure [15] are the results of the two ID runs. At 
this point, the 2D results are only slightly different from the ID results. However, in some region behind the shock, the 
gradient of entropy is negative. A negative entropy gradient region from ~ 70 km to ~ 100 km is also clearly visible in 
the entropy profile at 4.1 ms (left panel of Figure [TB)) . The white and black contour lines in the left panel of Figure [TBI 
denote s = 7.5 and 5 k baryon - , respectively. Convective instabilities can potentially develop in the region between 
~ 60 km and 90 km. In fact, convection does emerge quickly in our simulation. In the right panel of Figure [TBI we 
show the entropy profile at 9.8 ms after bounce. Also shown in the figure is the velocity field. High entropy plumes 
rise while some low entropy material falls down. Secondary Kelvin-Hclmholtz instabilities also develop in the shearing 
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Fig. 15. — Density (upper left), radial velocity (lower left), temperature (upper center), entropy per baryon (lower center), electron 
fraction (upper right), and net lepton fraction (lower right) as a function of enclosed mass at 1.5 ms after bounce. Here, radial velocity is 
the velocity in spherical radial direction. We show the results of the 2D cylindrical simulation (red solid lines), ID low-resolution spherical 
run (blue dashed lines), and ID high-resolution spherical run (green dotted lines). The radial profiles of the 2D results are generated via 
averaging. 
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Fig. 16. — Entropy at 4.1 (left) and 9.8 ms (right) after bounce for the 2D cylindrical core-collapse supernova simulation. The white and 
black contour lines in the left panel denote s = 7.5 and 5 k baryon -1 , respectively. The arrows in the right panel show the velocity field in 
the shocked region at 9.8 ms after bounce. It is shown that convection has developed in the shocked region at 9.8 ms after bounce. 

regions. The convection helps push the shock outwards because of the enhanced material energy transport from the 
inner hot region to the region behind the shock by the convective motion. Figure [17] shows the density, entropy, 
temperature, and electron fraction profiles at 20 ms after bounce. The shock has moved to ~ 200 km at 20 ms in the 
2D simulation (see also Figure [TS]). At 40 ms after bounce, the north pole and the south pole of the shock have moved 
to ~ 450km. In contrast, in the ID simulation in Section l4.2.1l (Figure [TU)) . the shock radius is ~ 85km at 20 ms and 
it never exceeds 141km. The comoving frame luminosities of v e , v e , and ty, neutrinos measured at 500 km for the 2D 
simulation are shown in Figure 1191 Also shown in Figure [19] are the results of the ID low-resolution run, which are 
nearly identical to those of the ID high- resolution results (Figure ITT]) . In the 2D results, there is also a large narrow 
pulse in L„ e that rises and drops at roughly the same rate. However, after ~ 7 ms past bounce, the 2D results are 
qualitatively different from the ID results because of the emergence of the convection. A major difference is that 
neutrinos are much less luminous in 2D than ID, because the temperature in 2D is cooler than that in ID due to the 
rapid expansion of 2D shock. 

Since this paper is not on core-collapse supernova explosions per se, the 2D simulation is stopped at ~ 55 ms after 
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Fig. 17. — Snapshot at 20 ms after bounce for the 2D cylindrical core-collapse supernova simulation. We show (a) logarithmic density, 
(b) entropy per baryon, (c) temperature, and (d) electron fraction, where the unit for density is g cm -3 . Only the inner 270 km of the 
model is shown here. 
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Fig. 18. — Shock radius vs. post-bounce time for the 2D core-collapse supernova simulation. We show the shock position in the northern 
hemisphere (red solid line), southern hemisphere (green dotted line), and equatorial plane (blue dashed lines) for the 2D simulation. Also 
shown is the result of the ID low-resolution simulation (black dash-dotted line). 




Fig. 19. — Comoving frame luminosities of v e , v e , and neutrinos measured at 500 km for the 2D cylindrical (red solid lines) and the 
low resolution ID spherical (black dashed lines) core-collapse supernova simulations. 

bounce. A detailed analysis of the 2D simulation is beyond the scope of this paper on algorithm, implementation, and 
testing. 

5. PARALLEL PERFORMANCE 

In paper II, we demonstrated that the gray radiation solver in CASTRO has very good scaling behavior on up to 
32768 cores via a weak scaling study. The multigroup radiation solver presented in this paper has similar good weak 
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Fig. 20. — Strong scaling behavior of 3D simulations on Hopper at NERSC. We show (a) wall clock time per coarse time step and (b) 
parallel efficiency in circles for five runs with the number of cores ranging from 1536 to 24576. The dotted line in panel (a) indicates perfect 
scaling. 

scaling behavior because the same module for solving Equation (|65j) is used in both solvers. Here we present a strong 
scaling study that was carried out on Hopper, a petascale Cray XE6 supercomputer at the National Energy Research 
Scientific Computing Center. We have performed a series of three-dimensional core-collapse supernova simulations 
in Cartesian coordinates with the same resolution on various numbers of cores. A hybrid MPI/OpenMP approach 
with 6 OpenMP threads per MPI task was used for these simulations. PFMG, a parallel multigrid solver from hypre 
(|Falgout fc Yandl200l \hvpre Code Projectl[20ll ). was used to solve the linear systems. In each run, we start the 3D 
simulation from the result of a ID simulation at about 10 ms before bounce. We use 12 energy groups for each species 
(v ei D e , and v^). The computational domain is a three-dimensional box with a size of 5000 km in each dimension. 
Three refinement levels (four total levels) are used with a refinement factor of 2 between level and level 1, and a 
refinement factor of 4 for other levels. The coarsest level has 256 3 cells. The finest level, which approximately covers 
the inner region of y/ x 2 + y 2 + z 2 < 165 km, has a cell size of 0.61km. Figure l20l shows the wall clock time per coarse 
time step for a series of runs with the number of cores ranging from 1536 to 24576. Also shown in Figure I2TJ1 is parallel 
efficient defined as 

*<M = fg, (in) 

where N is the number of cores, T/v is the wall clock time for the N-cove run, Nq is the number of cores in a baseline 
run (i.e., 1536-core run), and To is the wall clock time for the baseline run. In this strong scaling study, good scaling 
behavior is achieved on up to ~ 10000 cores. However, because this is a strong scaling study, the work load per 
core diminishes as more cores are used. For example, the average number of cells for each core in the 24576-core 
run is about 16 3 . An unfavorable consequence of small computational boxes is an increased communication overhead. 
Another unfavorable consequence is that more computation is wasted at domain boundaries due to larger surface area 
to volume ratios of smaller boxes. Hence, the degradation of parallel performance on more than 10000 cores is not 
surprising. 

6. SUMMARY 

In this paper, we have presented a multi-dimensional multigroup radiation hydrodynamics solver that is part of the 
CASTRO code. Block-structured AMR is utilized in CASTRO. In this paper, we focus attention on the single-level 
algorithms because the AMR algorithms have been presented in Papers I, II, and references therein. 

Our multigroup radiation hydrodynamics solver is based on a comoving frame formulation that is correct to order 
0(v/c), and uses a FLD approximation. Our mathematical analysis shows that the system we solve contains a 
hyperbolic subsystem that is the usual hydrodynamics system modified by radiation, a frequency space advection part 
that accounts for the Doppler shift of radiation due to the motion of matter, and a parabolic diffusion part. We have 
presented the mathematical characteristics of the hyperbolic subsystem. The eigenvalues and eigenvectors we have 
obtained are useful for characteristic-based Godunov schemes. We also described our treatment of the frequency space 
advection part. An implicit solver involving two nested iterations is presented. We have also presented two acceleration 
schemes that improve the convergence rate of the implicit solver. 

Extensive testing is presented demonstrating the accuracy and stability of CASTRO to solve multigroup radiation 
hydrodynamics problems. We have presented a scries of photon radiation test problems that cover a wide range of 
regimes. In addition to photon radiation problems, we have demonstrated the applicability of CASTRO in core-collapse 
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supernova simulations. 

Recently. iLentz et al.l (|2012[ ) have argued, via a series of ID simulations, that general relativistic effects and inelastic 
scattering of neut rinos on matter have si gnificant effects on the numerical modeling of core-collapse supernova explo- 
sions. In contrast. iNordhaus et al.1 (|2010f ) have concluded that spatial dimensionality has far more impact than general 
re lativistic effe c ts or microphysics such as inelastic scattering. Our ID simulations are consistent with the assessment 
of ILentz et al.l (120121) on the effects of inelastic scattering. However, our 2D simulation supports the assessment of 
INordhaus et al.l (|2010l) on the importance of spatial dimensionality. Our simulations have shown that multidimen- 
sional effects (e.g., convective motion in the region behind the shock and exterior to the nascent neutron star) appear 
less than 10 milliseconds after bounce, and 2D results are profoundly different from ID results. Thus, one should be 
cautious about assessment of various effects such as inelastic scatt ering that a r e bas ed on ID simulations that last 
several hundred milliseconds past bounce like the ones presented in ILentz et al.l (|2012l ). 

The main advantages of the CASTRO code are the efficiency due to the use of AMR and the accuracy due to the 
coupling of radiation force into the Riemann solver. Three-dimensional multigroup neutrino radiation hydrodynamics 
simulations of core-collapse supernovae with a good resolution using CASTRO can be carried out, if not now, in the 
near future with faster computers. In the strong scaling study with a total of 36 radiation groups and a cell size of 
0.6km on the finest level (Scction[5]), it took about 0.8 hours to advance one coarse time step (about 0.1 ms) on 12288 
cores of Hopper. Thus, it would take 800 hours on 12288 cores of Hopper to evolve to 100 ms. Admittedly, it is still 
very expensive to perform long-duration 3D simulations, but it is possible to study the first 100 ms after bounce and 
investigate the development of turbulent convection region via 3D simulations. 
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